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Abstract 

We derive the exact dual representation of the bosonic Hubbard 
model which takes the form of conserved current loops. The hardcore 
limit, which corresponds to the quantum spin-i Heisenberg antifer- 
romagnet, is also obtained. In this limit, the dual partition function 
takes a particularly simple form which is very amenable to numerical 
simulations. In addition to the usual quantities that we can measure 
(energy, density-density correlation function and superfluid density) 
we can with this new algorithm measure efficiently the order parame- 
ter correlation function, {aia}-), \i — j\ > 1. We demonstrate this with 
numerical tests in one dimension. 



1 Introduction 

The development of ever more efficient and sophisticated Quantum Monte 
Carlo (QMC) algorithms has greatly advanced our understanding of quan- 
tum many body systems and the various phases that exist in such models 
(superfluid, Mott insulators, Bose glasses, supersolids etc) and the transi- 
tions between them. Examples of such QMC algorithms include the "path 
integral algorithm" which was used to study the details of the Helium super- 
fluid transition in 2 and 3 dimensions Q, the World Line algorithm which 
was used to study the various phases of the bosonic Hubbard model. || 
Improvements of the World Line algorithm for hard core bosons (infinitely 



repulsive contact interaction) came in the form of cluster algorithms Q| where 
one updates many variables at a time. Such algorithms converge much faster 
and suffer much less from critical slowing down. The next improvement was 
to eliminate Trotter errors || which come from discretizing the imaginary 
time (i.e. inverse temperature) direction. These continuous imaginary time 
cluster algorithms H are the state of the art in hardcore boson simulations 
although their efficiency goes down in the presence of disorder and longer 
range interactions like next near neighbour. 

In the case of extreme soft-core high density bosonic models a differ- 
ent algorithm was developed based on duality tranformation. In this case, 
the amplitude of the order parameter is assumed constant, leaving only the 
phase, giving rise to a model of the XY variety called the Quantum Phase 
Model (QPM). The dual of this model is easily obtained when the Villain 
aproximation[|?J of the action is taken, leaving a model of interacting con- 
served integer current loops that can be quite easily simulated. § This loop 
algorithm was used productively by several groups to study the phases of 
the QPM with or without disorder. |8| 

One major disadvantage of the above algorithms is their inability to 
measure the correlation function of the order parameter, (ajaj) when \i— j\ 
is greater than 1. This quantity is very interesting since in phase transitions 
the behaviour of the order parameter is of prime importance. 

We will outline below how to perform the duality transformation ex- 
actly for the soft and hard core bosonic Hubbard models. Then, concentrat- 
ing on the hardcore case, we will construct the (exact) loop algorithm and 
demonstrate how it can be used to obtain the correlation function of the 
order parameter in addition to the the usual quantities of interest (energy, 
density-density correlation function, superffuid density). 

2 The Bosonic Hubbard Model 

We are interested in simulating models whose Hamiltonian has the form 

H = -t^(a\aj+a^ j a i )-^[iih i +Vo^ht+Vi^hih j +V2 ^ hih k . (1) 

{ij) i i (ij) {(ik)) 

In this equation, t is the transfer integral (the hopping parameter) which 

sets the energy scale (in the simulations we set it equal to 1), Vo, V\, and 
Vi are respectively the contact, the nearest neighbor, and the next near 
neighbor interactions, which we always take to be repulsive. On the square 
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lattice, the next near neighbor is chosen along the diagonal, while in one 
dimension the choice is the obvious one. dj and a] are destruction and 
creation operators on site i satisfying the usual softcore boson commutation 
laws, [o»,aj] = Sij, and hi = a\a i is the number operator on site i. In Eq. 
(1), (ij) and ((ik)) label sums over near and next near neighbors. The first 
term which describes the hopping of the bosons among near neighbor sites 
gives the kinetic energy. ^ = fx+Si is the site dependent chemical potential. 
For a clean system i.e. with no disorder, we take 5{ = 0. To include the 
effect of disorder, we take 5i to be a uniformly distributed random number 
— A < Si < A. We have chosen to disorder the system with a random 
site energy, which preferentially attracts or repels bosons to particular sites. 
Other choices are possible, such as disordering the hopping parameter, t, or 
some of the interaction strengths, Voi2- It is thought that these different 
ways of introducing disorder yield similar results. 

The quantum statistical mechanics of this system is given by the parti- 
tion function Z, 

Z = Tre- pB , (2) 

where (3 is the inverse temperature. The expectation value of a quantum 
operator is given by 

(O) = ±Tr(Oe-P B ). (3) 



3 Coherent State Representation 

We cannot do numerical simulations directly on this formulation of the par- 
tition function because it is in terms of operators. We must first find a 
c-number representation which can be incorporated into a simulation algo- 
rithm. There are several such representations. For example the wavefunction 
representation leads to what is known as the "path integral" method 0] and 
the occupation state representation leads to the World Line algorithm.^] 
Here we will use the coherent states representation, i.e. eigenstates of the 
destruction operator, 

a*|{*})=^(0K*}>. (4a) 
({d>}|al = (W|<m (46) 
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where the eigenvalues <p(i) are complex numbers that are defined on the 

sites, i, of the lattice. In terms of the more familiar occupation number rep- 
resentation vacuum (a\0) = 0,a^|0) = [1}) the coherent state |<3?) is defined 
as follows 

HZ}) = expfa-^^l + <K*k f )) 10). (5) 
With this normalization we obtain the inner product of two coherent states 

(mm) = ex P (^(i,^m) - \<t>*{m) - \rwm)\ w 

and the resolution of unity 



i=/n^i{*»«*}i- w 

i 

Now we write the partition function, Eq. (2), as 

Z = tr(e- 5H e- 5H e- 5H ...e- 5H ), (8) 



where 5 = (3/L T , L T is a large enough integer such that 5 <C 1. We express 
the partition function this way because we can now express the exponen- 
tials in Eq. (8) in a form suitable for easy (albeit approximate) evaluation. 
Between each pair of exponentials introduce the resolution of unity, Eq. 7. 
Using standard manipulations H we find 



where the action is given to first order in 5 by[10| 



S = J2 <f>*(r, r)A_ r 0(r, r) + 5 £ (r, r + 1), 0(r, r)}. (10) 

In Eq. 10, t denotes imaginary time, i.e. the inverse temperature direction 

(3. Also, H[(p*(r, r), (p(r, r — 1)] simply means that at the imaginary time r, 
we replace in the Hamiltonian, Eq. 1, the destruction operator a r by the 
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complex field <fi(r, r), and the creation operator, aj, by </>*(r, r + l).Q In this 
article our notation for the forward and backward finite difference operators 
are 

A M 0(r) = 0(r + £)-<Kr), (Ha) 

A_ M 0(r) = 0(r)-0(r-A), (116) 

It is well known|| that in the first term of Eq. 10 we must have A_ T and 
not A T although in the continuum limit one might argue that they both 
lead to the same result. Such arguments are incorrect. What is perhaps less 
strongly emphasized is that the following "approximation" , 

5 £ H[<f>*(r, r), 0(r, r - 1)] = H[^*(r, r), 0(r, r)] + 0(6/3). (12) 

r r 

which is often used in the literature, is not correct. The terms ignored are 
in fact of order f3, not 5/3, which in effect changes the Hamiltonian under 



consideration. [11] This can be easily illustrated using with the quantum 



harmonic oscillator. We will not use this approximation in what follows. 

4 The Duality Transformation 

In what follows we will take Vq = V\ = V% = in order to simplify the 
presentation. Nonzero values of these parameters can be dealt with in a 
straight-forward way. [11] Interestingly, the hardcore boson case (with no 
near and next near neighbor repulsion) is obtained directly from this case 
of zero interaction (see below). 

To effect the duality transformation on the partition function, Eqs. 9 
and 10, we follow the method of reference [12]. To this end we first write 
the coherent states field in terms of the amplitudes and phases: <j>(r, r) = 
a(r, T)e l9 ( r ' T \ Consequently, 

f [] d 2 cf>(r, t) - J J] a(r, r)da(r, r)d9(r, r), (13) 

^ r,r r,r 

where we are dropping irrelevant overall constant factors. The action 
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S = Y / U*(r,T)A_ T <t>(?,T) 

d \ 

-t<5^(^(f,r+l)0(r + A;,r)+^(r + fc,r + l)^(f,T)) , (14) 



becomes 



S = J2[ r)e- t9 ^ r) (a(r, T )e ie ^ T) - a(f, r - l)e ie ^ T - 1] ) 

d 

-t5j2 ( a (r, r + l)a(r + k, r y Wr+^)-0(r>+i)) 
k=i 

+ a(f, r)a(r + k, r + i) e i(e(^r)-0(r+fe,r+i))^ _ (15) 

Here we took the model to be in d space dimensions indicated by the index 

k = 1, d. k is & unit vector in the the kth. direction. 

We see in this equation that the phase, 6(f,T), which is a site variable, 
appears only as differences of near neighbor sites |l3|]. Thus we can see 
the beginings of an Xy-like model. The fact that we always have this 
combination of variables motivates us to change the variable of integration 
from the site variable 9{f, r) to the link or bond variables 0fc(r, r) = A/ t 9(f, r) 
and 6 t (t,t) = A T 0(f,r). In two space dimensions, the partition function 
thus becomes 




where e pi/p is the totally antisymmetric tensor in three dimensions. Note 
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that the product over plaquettes is over all space-space and space-time pla- 
quettes. The 5-functions are the only Jacobian of this variable change, fl~2|l 
Its geometrical interpretation is simple. Even though the model is now ex- 
pressed in terms of bonds, it is actually a site model. It is clear from the 
definitions of the link variables that if we sum them along any directed closed 
path the result is zero (mod 2tt). This is known as the Bianchi identity which 
is lost when the model is expressed in terms of bonds and needs to be en- 
forced as a constraint. The first set of ^-functions in the above equation 
enforces the "local" Bianchi identities, i.e. those due to topological!/ trivial 
loops. The second set enforces the "global" Bianchi identities, i.e. those 
due to topological!/ nontrivial loops in the imaginary time direction due to 
the periodic boundary conditions. These constraints have several interesting 
relationships to various geometrical aspects of the theory.^] Here we will 
simply exploit its relationship to the duality transformation. As was shown 
in Ref. [^] the dual variables are the Fourier conjugates to these constraints. 
In other words, the Fourier expansions of these 5- functions 
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i H fT e p( r »w A -" i f( ,? . T ) 



(17) 



(18) 



= E 

{n T (r)=-oo} 

immediately give the dual variables. In the case of the local Bianchi con- 
straints, Eq 17, the dual is the integer valued bond variable L(r, t), while 
in the global case, Eq 18, the dual variable is the integer valued field n T (r). 
Note that whereas l^(f, r) is a vector bond variable which depends on both 
coordinates r and r and which has components in the x, y, andr directions, 
the variable n T (r) is global and only points in the time direction. It depends 
only on the spatial coordinates and gives the value of the current flowing in 
the time direction from r = to r = L T where L T is the number of steps in 
the imaginary time direction. Jl2| Also note that the form in which the dual 



variable lp(f, t) appears is always e Atl/ pA_ I/ Z jU (r, r). So we define the local 
electric current 
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3p(r, t) = -e^pA^J^r, r). (19) 

Due to the totally antisymmetric tensor e, it is clear that the local integer 
current j p (f,r) is conserved. 

Substituting Eqs 17, 18 and 19 in Eq 16 we can integrate over the original 
variables, # M (r, r) and a(f,r). The details of the integration will be given 
elsewhere, [pi]] This leaves the partition function expressed only in terms of 
the dual variables n T , and s k . The new variable s k is positive semidefinite 
integer valued and is nothing but the dual to the amplitude field a. The 
partition function thus becomes: 

'Er 2 Sfc (r,r) / -i-r [n T (r) + j T {f,T)]l \ 



z= e (^^(n 

{™T,j M ,Sfc} 



n T {r) + jr(f,r) - M(f,r)]! 

• ■ ) ■ J 



J} h Mf, r)]l[s k (r, r) + Jk (f, r)]\ ) ' (2 ° } 



where 



M(r,r)= J2 ( s k{r,r)+s k (r-k,T)+j k (f-k,T). (21) 

k=l,2 

In this expression, which will be greatly simplified below, we allowed for the 

possibility of disorder in the chemical potential. Also, even though originally 
the values of n T ran from — oo to +oo (Eq. 18), the 6 T integrals impose the 
condition that n T is positive semidefinite. In fact we can easily prove|ll] 



that the total current traversing a bond in the time direction is nothing but 
the number of bosons traversing that bond. In addition, it is understood 
that all the arguments of factorials are positive semidefinite. This imposes 
several severe constraints on the allowed configurations which we will exploit 
to simplify the partition function. 

For simplicity, although this is not necessary, we will consider the no 
disorder case, i.e. fi(r) — > fi. In addition, we will consider the very important 
case of hardcore bosons where n T (r) + j T (r,T) can take only the values 
or 1. Combining this with the previously mentioned constraints on the 
arguments of the factorials, allows us to solve for the allowed electric loop 
configurations. Ill]] In this case the partition function simplifies drastically. 
For example for the one dimensional model it becomes 
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Z Q = E E e^^ n ^ x \t6)^ ]ix{x ' T)l , (22) 

{n T =0,l} {j M =0,±l} 

for the grand canonical ensemble and 

Z= (tS)^ lMx ' T)l , (23) 

{j M =0,±l} 

for the canonical case. These are the duals to the hardcore boson Hubbard 

model. The interpretation is very simple: The dual is a model of conserved 
integer current loops that take the values or 1 in the time direction and 
0, ±1 in the space directions. The partition function is a sum over all de- 
formations with each spatial hop (in this case x) costing t5. This is a very 
simple and appealing form and is very amenable to numerical simulations. 

The algorithm is now very simple. For example for the canonical case, we 
start with all local currents zero, and with the desired number of nonzero 
global time currents, n T , corresponding to the number of bosons in the 
system. Then we visit each plaquette and randomly choose to add to it 
a positive or negative elementary current loop (plaquette). This attempt 
is rejected if (i) it introduces negative time currents, or (ii) it introduces 
currents larger than 1. If this test is satisfied, we accept this current loop 
in accordance with detailed balance: If (t5) 2 is greater than or equal to a 
uniformly distributed random number we accept the change. 

5 Physical Quantities 

Now that we have the partition function and the algorithm we need expres- 
sions for the physical quantities we want to measure. The expression for the 
energy can be easily obtained from < E >= — ^lnZ. This gives 

(E) = -±(£\j x (x,T)\). (24) 

The numerical values for (E) obtained with the above algorithm and Eq 24 

are shown in Fig. 1 and compared with the results of exact diagonalization. 
We see that the agreement between the two results is excellent. 

The expression for the equal time density-density correlation function is 



9 



(n(x 1 )n(x 2 )} = ((rw(x 1 )+j T (xi,T))(n T (x2)+j T (x2,T))), (25) 

since the number of particles on site (x,r) is given by (n T (x) + j T (x,r)). 

We compare the numerical values of this correlation function with the exact 
values in Fig. 2. Again we see that agreement is excellent. 

The superfluid density is, in general, related to the winding number, W, 
of current configurations by|| 

^ fort 

ps = w (26) 

We have so far only discussed local moves in our Monte Carlo algorithm, 
and such moves can never change the winding number of a configuration. 
Therefore, the system is stuck in the winding number sector of the initial 
configuration. If the initial configuration has W = 0, this will give (W 2 ) = 
and therefore zero superfluid density. This obstacle can be overcome with a 
trick. Define the total x-current at a given time r by J x (t) = J2 X Jx( x > T ) an d 
calculate the Fourier transform of the current-current correlation function 

J(uj) = ^2^ WT (Uro)Uro + r)}. (27) 

t,to 

We can show|Q that J(uj — > 0) = W 2 which allows us to calculate the 

superfluid density p s . In Fig. 3 we show p s as a function (p c — p) where 
p c = 1 is the critical density at which the hardcore boson system becomes 
an incompressible Mott insulator and p is the boson density. From scaling 
arguments [|l4j] we expect p s ~ (1 — p) uz with vz = l. p5[ Our numerically 
obtained value on a small system is vz = 0.96 in excellent agreement with 
the theoretical values. 

These three quantities can be easily measured by existing QMC algo- 
rithms mentioned earlier. What is qualitatively new here is that we can 
measure the correlation function of the order parameter very easily. Let 
N(r) be the number of jumps of length r in a spatial direction (in this case 
x) in a given loop configuration. We can show[p"l|] 

-#-]lEl— X2\ 

1 



(a(x 1 )aHx 2 )) = ( (l + W(\xi-x 2 \ + l)), (28) 

P° 1=0 



where |xi — x 2 \ > 1- This quantity is shown in Fig. 4 for f3 = 0.5 and 4 for 
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L x = 8. We see that agreement with exact values is excellent for (3 = 0.5 
and gets worse for (3 = 4 as — X2 1 increases. The reason is that the exact 
diagonalization results sum over all winding number sectors whereas our 
simulation stays in the W = sector. At high temperature (/? = 0.5) the 
correlation length is short and the boundary conditions are less important 
so we get good agreement. At lower temperatures the correlation length 
is longer and consequently the boundary conditions become very important 
for small systems. For the more interesting larger system sizes the boundary 
conditions will become less important. 



6 Conclusions 

We have outlined how to perform the duality transformation exactly for 
the soft and hard core bosonic Hubbard models. In the hard core case the 
dual model is particularly simple. The dual partition function was used to 
construct a loop Quantum Monte Carlo algorithm and we showed that the 
numerical results agree very well with those of exact diagonalization for a 
system with 8 sites and 4 bosons. 

In particular we showed that with our algorithm we can calculate easily 
the order parameter correlation function, {a{x\)a) {xz)) which is very difficult 
to do with most other QMC algorithms. This interesting quantity is very 
important for the elaboration of the quantum phase transitions exhibited by 
the system. 

This representation of the Hubbard model has many common features 
with the World Line representation for which very efficient cluster algorithms 
have been constructed.^, [| It, therefore, seems possible to improve the effi- 
ciency of this algorithm in the same way. A cluster algorithm would improve 
the convergence properties and in addition would lift the restriction to zero 
winding number sector. [Q] This would greatly improve the calculation of the 
order parameter correlation function. This is curently under investigation. 

Note added. After the completion of this work we became aware of 
two other algorithms which allow to calculate the correlation function of 
the order parameters. These are the "worm" algorithm fi~6[| and the cluster 



algorithm with a new interpretation of the clusters. 1 17] 
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FIGURE CAPTIONS 

FIG. 1 The average energy per site as a function of (3. L x , Nf, are the number 
of lattice sites and bosons respectively. 

FIG. 2 The density-density correlation function as a function of distance. 
Open symbols are exact diagonalization values, full symbols are sim- 
ulation results. The lines are just to guide the eyes. 

FIG. 3 The superfluid density, p s , as a function of (1 — p). The theoretical 
value of uz is 1. 

FIG. 4 The correlation function of the order parameter, < a(x\)a^ (X2) > as 
a function of distance for two values of (5. 
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